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Abstract 


In order to model the thermal structure of polythermal ice sheets accurately, energy-conserving schemes and correct 
tracking of the cold-temperate transition surface (CTS) are necessary. We compare four different thermodynamics 
solvers in the ice sheet model SICOPOLIS. Two exist already, namely a two-layer polythermal scheme (POLY) and 
a single-phase cold-ice scheme (COLD), while the other two are newly-implemented, one-layer enthalpy schemes, 
namely a conventional scheme (ENTC) and a melting-CTS scheme (ENTM). The comparison uses scenarios of the 
EISMINT Phase 2 Simplified Geometry Experiments (Payne and others 2000 J. Glaciol. 46, 227-238). The POLY 
scheme is used as a reference against which the performance of the other schemes is tested. Both the COLD scheme 
and the ENTC scheme fail to produce a continuous temperature gradient across the CTS, which is explicitly enforced 
by the ENTM scheme. ENTM is more precise than ENTC for determining the position of the CTS, while the per¬ 
formance of both schemes is good for the temperature/water-content profiles in the entire ice column. Therefore, the 
one-layer enthalpy schemes ENTC and ENTM are viable, easier implementable alternatives to the POLY scheme with 
its need to handle two different numerical domains for cold and temperate ice. 
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1. Introduction 


Many glaciers and ice sheets are polythermal with disjoint cold and temperate domains, separated by the cold- 
temperate transition surface (CTS) (Blatter and Hutter 1991[ ). Both the Greenland and Antarctic ice sheets are 
Canadian-type polythermal, that is, they are mainly cold, except for distributed temperate layers at the base where 
strain heating is largest and where there is a geothermal contribution. It is thus important to model the thermody¬ 
namics of ice sheets correctly by distinguishing both domains and accounting for the transition conditions between 
them. 

Various methods allow one to model the thermodynamic conditions in ice sheets. Thus far, SICOPOLIS (Simu¬ 
lation code for POLythermal Ice Sheets; e.g., Greve |1997b| Sato and Greve |2012[ Greve and Herzfeldl 2013 URL 
www.sicopolis.net) is the only three-dimensional ice sheet model that employs the polythermal two-layer scheme. In 
this method, the temperature and water-content fields in the two domains, cold and temperate ice, are computed on 
separate numerical domains, and the transition conditions at the CTS are used to track its position. 

In most older ice sheet models (e.g., Huybrechts 1990[[^lov and Hutter| 1996t Payne and Dongelmans[ 1997 


Ritz and others| 1997| l, the cold-ice method was applied by resetting any computed temperatures that exceed the local 
pressure melting point to the local pressure melting point. While very simple, this means that energy is lost, and the 
water content in the temperate layer as well as the transition conditions at the CTS are ignored. The cold-ice method 
has, however, always been available in SICOPOLIS as an alternative to the polythermal two-layer method. 

Aschwanden and others (2012|l introduced a new, enthalpy-based approach for ice sheet thermodynamics. In this 


method, the thermodynamic fields of temperature in cold ice and water content in temperate ice are replaced by one 
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common thermodynamic field, enthalpjQ for both domains, and only one common field equation must be solved. 
However, the Stefan-type energy- and mass-flux matching conditions at the CTS, which are important for determining 
its position (|Greve|[r997a|l, are not included explicitly in the formulation of the enthalpy scheme by |Aschwanden and] 


others (20121. Following the terminology of Blatter and Greve (20151, we refer to it as the conventional one-layer 


enthalpy scheme. This scheme has already been used in a number of ice sheet and glacier models (Brinkerhoff and 


Johnson 2013 Golledge and others] 2013t Seroussi and others 2013 Wilson and Flowers 2013 Gilbert and others 

2014 K einer and others[ 20151. 

Two different conditions of the CTS must be distinguished. For melting conditions, cold ice flows across the 
CTS into the temperate layer, where water starts to accumulate due to strain heating along trajectories. The opposite 
situation, freezing conditions, occurs further downstream, where temperate ice flows across the CTS into the cold 
domain and the accumulated water content freezes out, releasing latent heat. For melting conditions, the temperature 
gradient and the water content are continuous across the CTS, while for freezing conditions, discontinuities of these 
quantities occur ( |Greve||T997a[ ). Since the CTS tends to be rather steep near the terminus, only a small area of the CTS 
is freezing, and therefore, in results of ice sheet models, freezing conditions usually only occur at very few isolated 
grid points (Greve 1997b[). 


Kleiner and others (20151 tested the implementation of the conventional enthalpy scheme for a Canadia n-type 
parallel-sided slab in one finite-difference and two finite-element ice sheet models (TIM-FD^, ISSM, COMice). Blat¬ 


ter and Greve (2015|l compared the performance of four different versions of the enthalpy scheme for a parallel-sided 


slab with a custom-designed finite-difference program. Besides the conventional enthalpy scheme, they considered 
the two-layer front-tracking enthalpy scheme (an enthalpy version of the poly thermal two-layer scheme mentioned 
above), the one-layer melting-CTS enthalpy scheme and the one-layer freezing-CTS enthalpy scheme. In the two latter 
schemes, explicit tracking of the melting or freezing CTS, based on the respective transition conditions at the CTS, 
has been added to the conventional enthalpy scheme. An important finding of these works was that the conventional 
one-layer enthalpy scheme can produce correct solutions for melting conditions at the CTS, provided that the numer¬ 
ical handling of the discontinuity of the enthalpy diffusivity across the CTS is done carefully. However, especially 
for finite-difference techniques. Blatter and Greve ( 2015| l concluded that it is safer to use the one-layer melting-CTS 
enthalpy scheme, which enforces the transition conditions explicitly. For freezing conditions, the conventional one- 
layer enthalpy scheme fails because it cannot handle the associated discontinuities of the thermodynamic fields, and it 
is thus imperative to enforce the transition conditions at the CTS explicitly, as it is done in the one-layer freezing-CTS 
enthalpy scheme. 

For this study, in addition to the previously existing polythermal two-layer and cold-ice schemes, we have imple¬ 
mented the conventional one-layer enthalpy scheme and the one-layer melting-CTS enthalpy scheme in the SICOPO- 
LIS model. For the reason given above, freezing conditions are not considered here. We attempt to test and verify 
these four schemes in SICOPOLIS, and in particular to test how the various schemes handle the melting CTS between 
cold and temperate ice for Canadian-type polythermal situations in ice sheets. Based on the results of [Blatter and| 
Greve (2015|l, we consider the polythermal two-layer scheme to be the most reliable method and thus use its results as 


benchmark solutions. In Sections |^and[^we give an overview of the theory of ice-sheet thermodynamics and describe 
the implementation of the various schemes in SICOPOLIS. Section]^ gives the set-up of the scenarios derived from 
the suite of EISMINT (European Ice Sheet Modeling INiTiative) Phase 2 Simplified Geometry Experiments (jPayne 


[and othei^ |2000j l used for this study. In Section]^ we discuss the results, focusing on the simulated thickness of the 
temperate ice layer. Sectionj^concludes the paper. 


2. Outline of ice-sheet thermodynamics 

2.1. Standard polythermal thermodynamics 


The standard description of the thermodynamics of polythermal ice masses, for which we follow largely Greve 


(1997a I, is based on the fields of absolute temperature T in cold ice and water content W in temperate ice. The 


evolution equation for temperature in cold ice is given by 


dT I d I dT\ Q 

— -I- V ■ grad r = — — \ k — + — , 
ot pc oz\ oz I pc 


( 1 ) 


* Owing to the incompressibility of ice, the enthalpy is identical to the internal energy. 
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where t denotes time, z the vertical spatial coordinate, v the three-dimensional velocity vector, p = 910 kg m“^ the 
ice density, k the temperature-dependent heat conductivity of cold ice and c the temperature-dependent heat capacity 
of cold ice. Also, Q - tr(t ■ D) is the volumetric strain heating, where t is the Cauchy stress tensor, D the strain-rate 
tensor, the middle dot (■) denotes tensor contraction and tr the trace of a tensor. Horizontal diffusion terms have been 


neglected, which can be justified by scaling arguments making use of the shallowness of ice sheets (e.g. Greve and 
|Blatter|[2009] l. 

Similar to Eq. Q. the evolution equation for water content in temperate ice reads 

dW . 15/ dW\ Q 

—+ v.g,adH'=--(v—) + - 


^ + vg,adr„J + --(«— 


( 2 ) 


where v is the water diffusivity in temperate ice (assumed to be constant) and L - 3.35 x 10^ Jkg^' the latent heat of 
fusion. The very small terms in the second line of the equation arise from the fact that the temperature in temperate 
ice is not constant, but equal to the melting temperature that depends on the local pressure p. 


Tmip) ^Tq-Pp, 


(3) 


where To = 273.15 K is the reference temperature and p - 9.8 x lO^^KPa ' the Clausius-Clapeyron constant for 
air-saturated glacier ice (Hooke 2005| l. As in Eq. ([T]|, horizontal diffusion terms have been neglected in the evolution 
equation (j^. 

As already mentioned in Sect.[^ for melting conditions at the CTS, the temperature gradient and the water content 
must be continuous across the CTS ( [Greve 1997a| l. If we mark values at the cold side of the CTS by plus (H-) 
superscripts, values at the temperate side by minus (-) superscripts, and denote the normal unit vector pointing into 
the cold side by n, this reads 

grad n- grad T' ■ n (4) 


and 


= 0 . 


(5) 


Eor freezing conditions, the situation is more complicated in that the temperature gradient and the water content 
are in general discontinuous across the CTS ( [Blatter and Hutter [1991 Greve 1997a| l. However, as already mentioned 
in Sect.[T] this situation usually occurs only on very small parts of the CTS, and therefore we do not consider freezing 
conditions in this study. 

At the ice surface, we prescribe the surface temperature as a Dirichlet-type boundary condition. At the ice base, 
three different cases must be distinguished. Eor a cold base, the geothermal heat flux determines the normal derivative 
of the temperature (Neumann-type boundary condition). Eor a temperate base with no temperate ice layer above, the 
temperature is equal to the pressure melting point. Eor a temperate base with a temperate ice layer above, in addition 
to that, a boundary condition for the basal water content is required (unless the water diffusivity v is neglected). In 
order to influence the solution as little as possible, we choose the Neumann-type zero-flux condition dWjdz - 0. 

In order to prevent the water content in temperate ice from rising to unreasonable levels due to the accumu¬ 
lated strain heating, a drainage model is required. We have kept the simple formulation previously implemented in 
SICOPOLIS, which assumes that any water exceeding the prescribed threshold of = 0.01 (= 1%) is drained 
instantaneously. The transport mechanism to the ice base remains unmodeled; however, the amount of drained water 
is added to the computed basal melt rate, and thus accounted for in the computation of the vertical velocity and the 
evolution of the ice thickness. 


2.2. Enthalpy method 

As an alternative to the poly thermal thermodynamics outlined in Sect. 2.1 Aschwanden and others (20121 intro¬ 
duced the enthalpy method. Its strength is that it combines the temperature and water content into a single thermody¬ 
namic field, the specific enthalpy, for which a single field equation holds. 
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The specific enthalpy h of ice at absolute temperature T and water content W is given by 


h(T, W) 


1 

f 


c(T)dT + LW. 


( 6 ) 


Inversely, the temperature and water content are not unique functions of the enthalpy because of the dependence of 
the melting temperature on the pressure (Eq. (0). We denote the inverse of Eq. for zero water content (W - 0) by 
T{h) and the enthalpy of ice at the melting point for zero water content by h^. 


TAP) 

h^(p)^h(T^(p),W^0)^ j c{f)df. 

To 


The temperature and water content can then be obtained by 


T(h,p) 


T{h), h<h^{p), 

Traip), h > h^ip) , 


W{h, p) 


0, h< h^ip ), 

(h-hUp))IL, h>h^(_p). 


(7) 


( 8 ) 

(9) 


The balance equation for enthalpy reads 


where 


-h V ■ grad h - 

ot 


_5 

dz 





h < h^ip), 


h > hj,p) 


( 10 ) 


( 11 ) 


are the enthalpy diffusivities in cold (subscript c) and temperate (subscript t) ice, respectively. Like in Eqs. ([T]i and 
horizontal diffusion terms have been neglected in Eq. ([T§, while vertical diffusion is present in both cold and 
temperate ice. 

The transition conditions for a melting CTS, Eqs. Q and Q, transform into physically equivalent enthalpy con¬ 
ditions. Inserting Eq. Q in Eq. (j^ and considering that the temperature on both sides of the CTS is equal to the 
pressure melting point Tu, yields the continuity of the enthalpy. 


( 12 ) 

On the cold side of the CTS, W - 0 holds everywhere, so that differentiating Eq. yields grad/j^ = c(Tj„)gTadT^. 
On the temperate side, differentiating Eq. 0 yields grad/Zj^j = c(7'ni) grad T^. Hence, the enthalpy form of Eq. 0 is 

grad ■ n - grad ■ n, (13) 


which expresses the continuity of the sensible heat flux. Note that the enthalpy gradient across the CTS is in general 
discontinuous (grad - n > grad h~ ■ n) because the enthalpy in the temperate ice layer becomes larger than away 
from the CTS due to the increasing water content. 

Eor freezing conditions, in general the enthalpy and the sensible heat flux are discontinuous (h^ < h^, grad -n < 
grad /z^ ■ n). The exact form of the transition conditions shall not be given here. 

By employing Eq. (j^, the boundary conditions for the ice surface and base, as well as the simple drainage model, 
given in the last two paragraphs of Sect. |2^ translate readily to the enthalpy formulation. However, for a temperate 
base with a temperate ice layer above, instead of the zero-flux condition for the water content dWjdz = 0 we use 
the zero flux condition for the enthalpy dhjdz - 0. Since the temperature in temperate ice is not constant (Eq. 0). 
these conditions are not equivalent, but the difference is small and both conditions are ad-hoc anyway, so that this is 
acceptable for the sake of simplicity. 
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3. Thermodynamics solvers in SICOPOLIS 


Previous versions of SICOPOLIS (3.1 and older) contained two different options for dealing with ice-sheet ther¬ 
modynamics. These are the polythermal two-layer scheme and the cold-ice scheme. Here, we added two more options 
based on the enthalpy method and the schemes described by Blatter and Greve| ( 2015| l, namely the conventional one- 
layer enthalpy scheme and the one-layer melting-CTS enthalpy scheme. Freezing conditions at the CTS are not 
considered. 


Blatter and Greve (20151 elaborated the schemes only for a one-dimensional problem (parallel-sided slab). Here, 


we extended them to three dimensions. SICOPOLIS uses terrain-following (“sigma”) coordinates, and the transforma¬ 
tion to these coordinates produces extra terms (e.g., Greve and Blatter] 2009 Sect. 5.7.1). However, since horizontal 
diffusion has been neglected in Eq. even in the transformed coordinates no mixing of horizontal and vertical 
derivatives of the enthalpy h occurs. This allows employing an implicit discretization scheme for the vertical deriva¬ 
tives and an explicit scheme for the horizontal derivatives in the same way as for the temperature and water-content 
equations ([T]| and (j^ ( [Greve |1997bl l. Thus, the numerical solution for each vertical profile is similar to the solution 
of the one-dimensional problem by an implicit scheme, and the horizontal advection terms of the three-dimensional 
equation play the role of additional, explicitly discretized source terms. 


3.1. Polythermal two-layer scheme (POLY) 

The polythermal two-layer scheme (scheme code; POLY) is the most sophisticated, but also the most complex 
method to simulate polythermal ice masses. It splits the computational domain numerically into two distinct layers of 
cold and temperate ice, and solves the evolution equations for temperature in cold ice (Eq. ([^) and water con tent in 
temperate ice (Eq. ([^) on two different grids ( [Blatter and Hutter| 1991| Greve 1997a[ Pettersson and others[ 2007| l. 
In this method, the CTS is hxed with the lower and upper boundaries of the cold and temperate domains, respectively, 
and thus can be tracked very precisely with the transition condition at the CTS. 

The method works for both melting and freezing conditions at the CTS; in particular, it can cope with the discon¬ 
tinuities of the temperature gradient and the water content that accompany freezing conditions. It is implemented in 
SICOPOLIS as an iterative trial-and-error procedure. Eor each time step and each ice column for which an internal 
CTS is detected, the temperature problem for the upper (cold-ice) region and the water-content problem for the lower 
(temperate-ice) region are solved repeatedly with different test positions of the CTS until all boundary and transition 


conditions are fulhlled. This is discussed in more detail by Greve (1997b last paragraph of section 2). Eor the simula¬ 
tions discussed in this study, we disregard freezing conditions and apply the POLY scheme only for assumed melting 
conditions at the CTS. The ability to turn off the freezing conditions is a newly-implemented option for the POLY 
scheme. 


3.2. Cold-ice scheme (COLD) 

The cold-ice scheme (scheme code: COLD) is the most simple, yet physically incorrect way to deal with poly¬ 
thermal conditions in ice sheets. SICOPOLIS applies the COLD scheme by solving the temperature equation ([T]) for 
the entire domain and resetting any temperatures T exceeding the local pressure melting point (Eq. (|^) to T -T^. 
The excess heat of the temperate layer becomes lost, so that the scheme does not conserve energy. 


3.3. Conventional one-layer enthalpy scheme (ENTC) 

The conventional one-layer enthalpy scheme (scheme code; ENTC) corresponds to the enthalpy method by [As-[ 


chwanden and others (2012[l. The enthalpy equation ([T0| is solved for the entire polythermal domain on a single 


numerical grid, and the CTS is diagnostically determined as the uppermost grid point of the temperate layer, that is, 
the uppermost grid point for which h > h^n holds. The transition conditions at the CTS, Eqs. ( [T^ and ( [T3] l for melting 
conditions or their more complicated counterparts for freezing conditions, are not enforced explicitly. 
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3.4. One-layer melting-CTS enthalpy scheme (ENTM) 


The one-layer melting-CTS enthalpy scheme (scheme code: ENTM) by [Blatter and GreTO ( 2015| l differs from the 
conventional one-layer enthalpy (ENTC) scheme in that the continuity condition (|13|l at a melting CTS is enforced 
explicitly. This is achieved for each vertical column possessing an internal CTS (that is, a non-vanishing thickness 
of temperate ice at the base) in two steps. In a predictor step, the enthalpy profile is computed for the entire vertical 
column like in the ENTC scheme. The CTS is then positioned at the uppermost point in the temperate layer. In 
a corrector step, the enthalpy profile in the upper, cold layer only is recomputed by using condition ( [T3| ) as lower 
boundary condition. The updated enthalpy profile then consists of the predictor step for the temperate layer and the 
corrector step for the cold layer. Eor the technical details of the scheme see Blatter and Greve ( 2015j l. 

The implementation of both the ENTC and the ENTM scheme in SICOPOLIS has the same topological restriction 
as the POLY scheme. It only allows for Canadian-type polythermal conditions, that is, only one CTS can exist in each 
column of ice, and the cold layer is always situated above the temperate layer. However, this is not a fundamental 
restriction for the enthalpy schemes. In principle, they also allow dealing with more complex and/or changing topolo¬ 
gies such as temperate ice at the surface, prescribed water content instead of prescribed temperature as a boundary 
condition at the surface, bubbles of temperate ice within cold ice etc. 


4. Experiment set-up 


We tested the four different schemes describe d in Sect, [^with m odified versions of the EISMINT Phase 2 Simpli¬ 
fied Geometry Experiments A and H as defined in Payne and others ( 2000) l. The computational domain in the horizon¬ 
tal plane is a square of 1500 by 1500 km, spanned by the Cartesian coordinates .r = 0... 1500 km, y = 0... 1500 km. 


The numerical values of the dynamical and thermodynamical parameters are given in Payne and others (20001 (see 


also Sect.j^. The experiments were performed using SICOPOLIS version 3.2-dev (revision 471) with three different 
combinations of horizontal grid resolution Ax (- Ay) and time step At, namely (Ax, Af) = (10 km, 2 a), (5 km, 2 a) 
and (10km, 20 a). Eor the 10-km resolution, this leads to 151 x 151 grid points (indices 1 = 0... 150, j - 0... 150) 
in the horizontal domain, while for the 5-km resolution the domain is covered by 301 x 301 grid points. The stan¬ 
dard vertical resolution is 81 grid points in the upper domain (terrain-following vertical coordinate fc = 0... 1, index 
kc = 0... 80) and, for the POLY scheme, 11 grid points in the lower domain (terrain-following vertical coordinate 
= 0... 1, index kt = 0... 10). Eor the other three schemes, the lower domain is not used. Eor the (10km, 2a) 
combination, we also consider a five times higher vertical resolution (abbreviated as “hvr”) with 401 grid points in the 
upper and 51 grid points in the lower domain. 

Both the heat capacity c and the heat conductivity k depend significantly on temperature (e.g. Greve and Blatterj 
2009) . The implementation of all thermodynamics schemes in SICOPOLIS accounts for that possibility by allowing 
both quantities to vary in space and time. This is particularly relevant for discretizing the vertical diffusion terms 
in Eqs. (0, (ig. and ( [T0| that contain k within d/dz terms. Howeve r, in the EISMINT set-up used here the values 
are taken as constants, c = 2009 Jkg *K ' and x = 2.1 Wm 'K ' (Payne and others 2000 1 . The diffusive flux in 
the temperate layer is likely very small for small water contents; however, virtually no information on the value of 
the water diffusivity v is available. In SICOPOLIS, a value of v = 10 ®kgm * s ' is implemented for reasons of 
numerical stability rather than for physical reasons. This leads to kt = v/p a; 1.1 x 10“^ m^ s“*, while kc = x/(pc) =» 
1.1 X 10“^ m^ s“', thus a three orders of magnitude difference. 


4.1. EISMINT experiment A1 

The experiment is started with no-ice conditions and runs for 200 ka to reach a steady state. The ice sheet rests on 
a flat, rigid base (no isostatic adjustment). The boundary conditions are circularly symmetric with respect to the center 
of the domain. The prescribed surface temperature increases linearly with distance from the center, with a minimum 
value of 238.15 K at the center (summit) and a horizontal gradient of 1.67 x 10“^ Kkm^'. Similarly, the surface mass 
balance is prescribed by a piece-wise linear function with a maximum accumulation rate of 0.5 ma * in the interior 
of the ice sheet (from the center to 400 km away from the center), and a horizontal gradient of -10“^ ma ' km ' 
beyond that, which leads to an equilibrium line 450 km from the center and an ablation zone further out. With these 
assumptions, the evolving ice sheet thickness does not feed back on the surface mass balance. At the ice base, no-slip 
conditions are assumed everywhere, and the geothermal heat flux is 42 mW m“^. 
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There are two differences compared to the original EISMINT experiment A. The first one is the higher horizontal 
resolution; we use 10 and 5 km (see above) compared to the original 25 km. The second one is that, following 
|Lliboutry and Du\^ ( |1985| l, we employ a water-content-dependent rate factor At for temperate ice, 

At(VT)=Ao=cX(l-Hl.81251T[%]), (14) 

where Aq^c = 4.529 x 10“^^ s“^ Pa“^ is the rate factor for 0°C (relative to pressure melting) that results from the 
original EISMINT set-up. 


4.2. EISMINT experiment A2 

The set-up of experiment A2 is essentially the same as that of experiment Al; however, the water-content- 
dependent rate factor ( [l4| ) is not used. Instead, the original EISMINT set-up is used, which only considers the 
dependence of the rate factor on temperature. Physically, this means that the lubrication effect (through enhanced 
shear) of the near-basal temperate ice due to its water content is ignored. 


4.3. EISMINT experiment HI 

Experiment HI replaces the no-slip condition of Al by basal sliding. Otherwise, the set-up is the same. The 
sliding law is most simple and relates the sliding velocity linearly to the basal drag. In the original experiment H, the 
sliding parameter has the value 10“^ ma ' Pa"'. However, this does not produce any basal temperate layers except 
for the COLD scheme due to the strong, downward advection of cold surface ice. Since we want to focus on the 
temperate ice, we decided to use a modified set-up, with the sliding parameter reduced by a factor 10 to the value 
Cl = lO-^'ma-'Pa"'. 

Eurther, in the EISMINT set-up, basal sliding is assumed to occur wherever the ice base is at the pressure melting 
point, while for a cold base no-slip conditions prevail. However, when used with the shallow ice approximation 
(that is employed by SICOPOLIS), this binary switch is known to produce a singularity of the vertical velocity field 
that prohibits proper convergence of the numerical solution ( |Bueler and Brown 20091. In order to circumvent this 
pr oblem, we regularize the transitio n by allowing for sub-melt sliding in the form by Greve (20051 (originally proposed 
by Hindmarsh and Le Meur (2001The sliding parameter Cb depends on the basal temperature T^ (in °C, relative 
to pressure melting) via 

Cb = C"exp(r'/y), (15) 

where y - 1°C is the sub-melt-sliding parameter. Equation ( [T5] l has the effect that, rather than stopping abruptly, basal 
sliding decays exponentially as the basal temperature falls below the pressure melting point. 


5. Results 

Eor experiment Al, run with the polythermal two-layer (POLY) scheme and the parameters Ax = 10km, Af = 2 a, 
Eig.[2depicts the simulated steady-state ice thickness (= surface topography) at the end of the simulation (f = 200 ka). 
The ice sheet reaches a volume of 2.109x10® km^, covers an area of 1.046x10® km^ and reaches a maximum thickness 
of 3.687 km. The melt fraction (fraction of basal ice at the pressure melting point) is 67.5%. 

Eor the same set-up, Eig.|^ shows four steady-state temperature profiles (at t - 200 ka) in the bottom-most 80 m 
of the ice. They are all located near the margin, ~ 520km away from the center, where the largest thicknesses of the 
temperate layer occur. The positions of the profiles are indicated in Eig. The profiles show the typical patterns 
of the four different schemes concerning the position of the CTS, how they match the transition conditions at the 
CTS and the influence of the discrete grid. Evidently, the POLY scheme produces the best solutions. The computed 
temperature profiles are both continuous and smooth across the CTS, thus fulfilling the transitions conditions for 
melting conditions. This justifies our decision (already stated in the last paragraph of the introduction) to use the 
solutions computed with the POLY scheme as a reference to assess the performance of the other three schemes. 


^An alternative way of regularizing the singularity at the spatial onset of hasal sliding is to replace the shallow-ice stress balance hy a higher- 
order one including membrane stresses, such as the first order approximation (Blatterj |1995) or a shallow-ice-shelfy-stream hybrid (Bueler and| 
|Brown||2Q09[|Pollard and DeConto||2012|. 
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Figure 1: EISMINT experiment Al: Simulated ice thickness (= surface topography) after 200ka model time. Polythermal two-layer (POLY) 
scheme, grid resolution Ajc = 10 km, time step At = 2 a. The positions of the four profiles shown in Fig.[^are maiLed by their respective grid 
indices Within grid resolution, they ai‘e at the same distance from the center [(35,42), (42,35) and (115,108) 518.6km away, (112,112) 

523.3 km away]. 


In the depicted four profiles, the COLD scheme overestimates the thickness of the temperate layer by a factor 
-'3.5. This overestimation, which was already reported by Greve (1997b I for simulations of the Greenland ice sheet 
with the SICOPOLIS model, is paralleled by a violation of the transition condition 0: the temperature gradient 
across the CTS is not continuous (most pronounced in Fig.|^). 

The profiles obtained with the ENTC scheme also overestimate the thickness of the temperate layer, but to a lesser 
extent than the cold-ice method (factor ~ 1.5-2). Like for the COLD scheme, the transition condition 0 (or its 
equivalent enthalpy form, Lq. ( [T3] l) is not fulfilled at the CTS. 

The temperate layer thicknesses found by the LNTM scheme are most similar to those produced by the POLY 
scheme. Lurther, the temperature profiles all show a continuous gradient across the CTS, that is, they fulfill the 
transition condition (|4 1 (or the equivalent form ([T3]l). However, the price to pay for its enforcement by the cotTector 
step described in Sect. 3.4| is that the temperature itself shows a slight discontinuity of ~ 0.02-0.03 K. Ultimately, this 
results from the fact that, in a one-layer scheme, the positioning of the CTS is naturally limited by the grid resolution 
(while, in the two-layer POLY scheme, arbitrary precision can be achieved). 

Comparing the temperature profiles in the cold-ice layer above the CTS shows that, despite the differences in the 
position of the CTS, the results produced by the POLY, LNTM and LNTC schemes are very close to each other. In 
detail, for all four profiles, the LNTM scheme is slightly warmer, while the LNTC scheme is slightly colder than the 
POLY scheme. By contrast, the COLD scheme produces notably (~0.4K) higher temperatures. 

Lor the POLY, LNTM and LNTC schemes, the majority of grid points within the temperate layer reach a water 
content of VT = VTmax = 0.01, which is the prescribed threshold value of the simple drainage model (see the last 
paragraph of Sect. 2.11. Therefore, we refrain from showing the water content in a separate figure. The COLD 
scheme is an exception because it does not account for any water content; in other words, the water content is zero 
everywhere. 

In the following, we focus on the volume and thickness distribution of the simulated temperate ice layers. This is 
because (a) the thermal conditions at and near the ice base (where most of the shearing takes place) are most relevant 
for ice flow, and (b) within the temperate ice, the water content usually reaches the value VTniax = 0.01 (see above). 
Ligure shows the evolution of the temperate ice volume for experiment Al, run with the four thermodynamics 
schemes, the grid spacing Ax - 10 km, the two time steps Af = 20 a and 2 a, and both the standard and high vertical 
resolution (“hvr”, see Sect. set-ups. The period is limited to the first 20ka of the total 200ka model time, during 
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Figure 2: EISMINT experiment Al: Vertical temperature profiles after 200ka model time at four different positions near the diagonal symmetry 
axis of the ice sheet (marked by their respective grid indices (/, j), see Fig. ^ in the zone of maximum thickness of the temperate layer. Grid 
resolution Ax = 10 km, time step Af = 2 a. Black dashed: POLY scheme, gray dashed: COLD scheme, gray solid: ENTC scheme, black solid: 
ENTM scheme (see Sect. I^for the scheme codes). The horizontal lines mark the CTS positions. 


which most of the changes take place. 

For all shown cases, the temperate ice volume starts forming 6-7 ka after the initiation of the ice-sheet build-up, 
goes through a weak maximum and then converges towards a steady-state value. The main difference between the 
larger (20 a, panel a) and the smaller (2 a, panel b) time step is the behavior of the COLD scheme. Compared to the 
POLY scheme, for Af = 20 a it overpredicts the temperate ice volume by more than a factor 6, while for Af = 2 a 
the overprediction factor is ~ 4. The results of the POLY, ENTM and ENTC schemes do not change much; ENTC 
produces ~2x more and ENTM ~ L4x more temperate ice than POLY. Employing the high vertical resolution (panel 
c) reduces the discrepancy between the POLY and ENTM schemes to less than 10%, and the overprediction by the 
ENTC scheme to a factor ~ 1.4, while there is little influence on the poor performance of the COLD scheme. 

Eor both time steps and the standard vertical resolution (panels a, b), the ENTC and ENTM schemes show some 
high-frequency noise in the evolution of the temperate ice volume (amplitudes less than 2% of the total temperate 
ice volume), while the results produced by the POLY scheme and the COLD scheme are smoother. The noisiness is 
clearly reduced (but still visible) for the high vertical resolution (panel c). This makes the one-layer approach of ENTC 
and ENTM, with its positioning of the CTS at a discrete grid point (uppermost grid point of the temperate layer, see 
Sects. 3.3 and 3.41 a likely cause. However, the one-dimensional simulations carried out by |B latter and Grev^p015| 
for a parallel-sided slab did not show any noise, so that its occurrence it also related to the three-dimensionality of the 
problem considered here (see below). 

Eigures^l^show maps of the steady-state thickness of the temperate layer for experiment Al, all four thermo¬ 
dynamics solvers and the same set-ups as in Eig. The decreasing amount of temperate ice in the order COLD > 
ENTC > ENTM > POLY is clearly visible. The larger time step in Pig. [^compared to Pig. [^mainly affects the results 
of the POLY and COLD schemes. For the COLD scheme, this was already discussed above. For the POLY scheme. 
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Figure 3: EISMINT experiment Al: Evolution of the temperate ice volume over the first 20ka model time. Black dashed: POLY scheme, gray 
dashed: COLD scheme, gray solid: ENTC scheme, black solid: ENTM scheme (see Sect. I^for the scheme codes). Grid resolution A;c = 10 km. a: 
Time step At = 20a. b: Time step At = 2 a. c: Time step At = 2 a, high vertical resolution (see Sect. |4]for details). 


a numerical instability occurs in the form of a waviness of the thickness of the temperate layer oriented diagonally to 
the X- and y-axes. The smaller time step used in Fig. [^resolves this problem. 

The larger temporal noise produced by the two enthalpy schemes (as discussed above) finds its counterpart in 
larger spatial noise (compare the respective panels c and d to a and b). The two types of noise are strongly related 
to each other. Both are largely unaffected by the time step (Fig. |5]vs. Fig. 1^1, but significantly reduced by the higher 
vertical resolution (Fig.|^vs. Fig.|^. For the small-time-step/high-vertical-resolution set-up used in Fig.|^ the results 
of the POLY and ENTM schemes are very similar to each other. 

Further, in all cases, the thickness of the temperate layer tends to increase gradually with distance from the center 
of the ice sheet, and then decreases sharply near the ice margin. This confirms for our new simulations what was 
already claimed in Sect. and it is the reason for the fact that freezing conditions at the CTS can hardly be resolved, 
thus justifying that we ignored them for this study. Owing to the incompatible symmetries of the ice sheet geometry 
(circularly symmetric) and the numerical grid (aligned to the x- and y-axes), the thin rings of largest temperate layer 
thicknesses near the margin are in all cases somewhat discontinuous. 

A synopsis of the main results of the 16 set-ups of experiment Al (all set-ups shown in Figs. |4||^ plus the high- 
horizontal-resolution case Ajc = 5 km. At - 2 a) is given in Table [T] The volume of temperate ice is in the range 
of ~ 0.25-1% of the total ice volume and, as discussed above, varies strongly across the different thermodynamics 
schemes. By contrast, the melt fraction (area ratio of basal temperate ice to total ice) is much larger (~2/3) and varies 
very little across all 16 set-ups. For the total ice volume and area, the strongest systematic (but still small) influence 
arises from the horizontal resolution; On average. Ax = 5 km produces an ice sheet with 0.7% more volume and 0.5% 
more area than Ax = 10 km. By contrast, the influences of the thermodynamics solver, the time step and the vertical 
resolution are unsystematic and even smaller. 

In Table 1^ the computing times of the four thermodynamics schemes are compared for Experiment Al with 
Ax = 10 km and At = 2 a. The differences between the schemes are significant. As expected, the most simple, 
but physically inadequate COLD scheme is by far the fastest, with a gain of more than 40% compared to the POLY 
scheme. The differences between the POLY, ENTM and ENTC schemes are smaller. ENTC is about 10% faster than 
POLY, while ENTM is about 6% slower. The reason for the relative slowness of the enthalpy schemes, despite greater 
simplicity compared to the POLY scheme, is that, in the SICOPOLIS implementation, in each time step the computed 
enthalpy field is converted back to temperature and water content. This must be done three-dimensionally and is thus 
costly. Codes that are designed as enthalpy-based from the outset (whereas the thermodynamics of SICOPOLIS is 
natively based on temperature and water content) can probably avoid these conversions, which opens the possibility 
for greater efficiency of the enthalpy schemes. 
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Figure 4: EISMINT experiment Al: Thickness of the basal temperate ice layer after 200 ka model time, a: POLY scheme, b: COLD scheme, c: 
ENTC scheme, d: ENTM scheme (see Sect.[^for the scheme codes). Grid resolution Ax = 10 km, time step At = 20 a. 
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Figure 5: Same as Fig.[^ but grid resolution Ax = 10 km, time step Af = 2 a. 
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Figure 6: Same as Fig.[^ but grid resolution Ajc = 10 km, time step At = 


2 a, high vertical resolution (see Sect. I^for details). 


13 


















Exp. Al set-up 

Volume 

Area 

Temp, volume 

Melt 

(scheme, Ajc, Af) 

(10^ km^) 

(10® km^) 

(lO^km^) 

fraction 

POLY, 10 km, 20 a 

2.105 

1.045 

0.537 

0.675 

COLD, 10 km, 20 a 

2.109 

1.043 

3.319 

0.671 

ENTC, 10 km, 20 a 

2.097 

1.044 

0.968 

0.673 

ENTM, 10 km, 20 a 

2.113 

1.046 

0.775 

0.674 

POLY, 10 km, 2 a 

2.109 

1.046 

0.501 

0.675 

COLD, 10 km, 2 a 

2.120 

1.044 

2.076 

0.675 

ENTC, 10 km, 2 a 

2.096 

1.044 

1.008 

0.676 

ENTM, 10 km, 2 a 

2.110 

1.046 

0.718 

0.674 

POLY, 5 km, 2 a 

2.123 

1.050 

0.466 

0.675 

COLD, 5 km, 2 a 

2.133 

1.049 

2.009 

0.674 

ENTC, 5 km, 2 a 

2.110 

1.050 

0.977 

0.675 

ENTM, 5 km, 2 a 

2.124 

1.050 

0.689 

0.675 

POLY, 10 km, 2 a (hvr) 

2.105 

1.045 

0.559 

0.695 

COLD, 10 km, 2 a (hvr) 

2.120 

1.044 

2.159 

0.696 

ENTC, 10 km, 2 a (hvr) 

2.101 

1.044 

0.790 

0.694 

ENTM, 10 km, 2 a (hvr) 

2.105 

1.045 

0.595 

0.694 


Table 1: EISMINT experiment Al: Mean values between t = 190ka and 200 ka for the ice volume, the ice area, the volume of temperate ice and 
the melt fraction (fraction of basal ice at the pressure melting point). See Sect. [^for the scheme codes. The abbreviation “hvr” means “high vertical 
resolution” (see Sect. I^for details). 


Exp. Al set-up 

Computing 

(scheme. Ax, Af) 

time (hrs) 

POLY, 10 km, 2 a 

7.7 

COLD, 10 km, 2 a 

4.2 

ENTC, 10 km, 2 a 

7.0 

ENTM, 10 km, 2 a 

8.2 


Table 2: Computing times for EISMINT experiment Al with the reference set-up {Ax = 10 km. At = 2 a), run with the Intel Fortran Compiler 15.0.3 
for Linux (optimization options -xHOST -03 -no-prec-div) on a 12-Core Intel Xeon E5-2697 v2 (2.7 GHz) PC under openSUSE 13.1 (64 bit). See 
Sect.[3]for the scheme codes. 


14 











Thickness of temperate layer (m) 

0 20 40 60 80 100 


750 


_550 

E 

^350 


150 

750 


550 

"e 

^350 


150 

150 350 550 750 150 350 550 750 

X (km) X (km) 



Figure 7: Same as Fig.[^ but for EISMINT experiment A2, grid resolution Ax = 10 km, time step Af = 2 a. 
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Figure 8: Same as Fig.[^ but for EISMINT experiment HI, grid resolution Ax = 10 km, time step Af = 2 a. 
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Exp. Al set-up 

Volume 

Area 

Temp, volume 

Melt 

(scheme. Ax, Af) 

(10® km^) 

(10® km^) 

(lO'^km^) 

fraction 

POLY, 10 km, 2 a 

2.109 

1.046 

0.501 

0.675 

COLD, 10 km, 2 a 

2.120 

1.044 

2.076 

0.675 

ENTC, 10 km, 2 a 

2.096 

1.044 

1.008 

0.676 

ENTM, 10 km, 2 a 

2.110 

1.046 

0.718 

0.674 

Exp. A2 set-up 

Volume 

Area 

Temp, volume 

Melt 

(scheme. Ax, Af) 

(10® km^) 

(10® km^) 

(lO^km^) 

fraction 

POLY, 10 km, 2 a 

2.124 

1.044 

1.245 

0.675 

COLD, 10 km, 2 a 

2.120 

1.044 

2.076 

0.675 

ENTC, 10 km, 2 a 

2.120 

1.044 

1.956 

0.675 

ENTM, 10 km, 2 a 

2.123 

1.044 

1.273 

0.676 

Exp. HI set-up 

Volume 

Area 

Temp, volume 

Melt 

(scheme. Ax, Af) 

(10® km^) 

(10® km^) 

(lO'^km^) 

fraction 

POLY, 10 km, 2 a 

2.094 

1.045 

0.079 

0.668 

COLD, 10 km, 2 a 

2.090 

1.041 

1.149 

0.659 

ENTC, 10 km, 2 a 

2.087 

1.045 

0.365 

0.664 

ENTM, 10 km, 2 a 

2.096 

1.045 

0.176 

0.664 


Table 3: EISMINT experiments Al, A2 and HI: Mean values between t = 190ka and 200ka for the ice volume, the ice area, the volume of 
temperate ice and the melt fraction (fraction of basal ice at the pressure melting point). See Sect. I^for the scheme codes. 


For experiment A2, the water-content-dependent rate factor for temperate ice has been replaced by a constant rate 
factor like in the original EISMINT set-up (Sect. [4~2l l. The resulting steady-state thickness of the temperate layer for 
all four thermodynamics solvers and the combination Ax = 10 km, Af = 2 a is shown in Fig.|^ and an overview of the 
main results is given in TableThe COLD scheme is unaffected by the difference between experiments Al and A2 
because it does not account for any water content, and consequently the results are the same. By contrast, the POLY, 
ENTM and ENTC schemes produce roughly two times thicker temperate ice layers for experiment A2. The reason 
is that the simulated temperate ice is stiffer, which reduces the advective transport of cold ice towards the base, so 
that the growth of the temperate ice layers is favored compared to experiment Al. It is interesting to note that, for 
experiment A2, the temperate ice volumes computed by the POLY and ENTM schemes agree very well (within 3%) 
for the standard vertical resolution employed here. For experiment Al, a similarly close agreement is only achieved 
with the high vertical resolution (see discussion above). 

There is no notable influence of experiment A2 on the melt fraction. The total ice area is essentially unaffected as 
well, while the total ice volume is systematically slightly (on average ~0.8%) larger than in experiment Al. 

We discussed above that, for experiment Al, the POLY scheme with the large-time-step setting Ax = 10km, Af = 
20 a produces a temperate-layer thickness that shows a waviness oriented diagonally to the x- and y-axes (Fig. |^). 
For experiment A2 and Ax = 10 km, Af = 20 a, this instability does not occur, and the thickness of the temperate 
layer (not shown) is very similar to that obtained for the setting Ax = 10km, Af = 2a shown in Fig. |^. Hence, 
the instability is linked to the enhanced fluidity in the temperate layer due to the water-content-dependent rate factor 
employed in experiment Al. 

For experiment HI, basal sliding over temperate- and near-temperate-based areas has been added to the set-up of 
experiment Al (Sect. [4^. Like for experiment A2, we only discuss the combination Ax = 10km, Af = 2 a. Figure]^ 
depicts the steady-state thickness of the temperate layer for the four thermodynamics solvers, and Table gives an 
overview of the results. Comparing Figs. and as well as the entries in Table [^reveals that the addition of some 
basal sliding reduces the volume of temperate ice strongly. This is a consequence of increased advection of cold 
surface ice downward and outward, and of reduced strain heating. Like for experiments Al and A2, the computed 
temperate ice volume decreases in the order COLD > ENTC > ENTM > POLY. However, even the ENTM scheme 
overpredicts the temperate ice volume by more than a factor 2 compared to the POLY scheme. This is because the 
thicknesses are so small that any one-layer scheme does not resolve them well (unless a very high vertical resolution 
is employed), whereas the performance of the two-layer POLY scheme is not limited by this problem. 
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Figure shows further that in the distribution of temperate ice obtained with the COLD scheme some spokes 
emerge. The spokes are also present in similar form in the results of the three other schemes when the basal tempera¬ 
ture is plotted instead of the thickness of the temperate ice layer (not shown). This phenomenon was already observed 
in the original article on the EISMINT Phase 2 Simplified Geometry Experiments ( Payne and others] 20001. It is most 
likely a consequence of the combined effects of a thermoviscous flow instability and the symmetry of the numerical 
grid (Bueler and others 2007) 1. 

Since the sliding coefficient in experiment HI is small, the melt fraction and the total ice volume only decrease 
slightly compared to those of experiment Al, and the ice area is virtually unchanged. This underlines the very high 
sensitivity of the temperate ice layer to changed flow conditions. For the original experiment H with a 10 times larger 
sliding coefficient (results not shown), the large-scale response of the ice sheet is more pronounced. For this scenario, 
only the COLD scheme produces a temperate layer of non-zero thickness, while the other schemes merely produce a 
temperate base. 


6. Discussion and conclusion 

We tested four different thermodynamics solvers in the ice sheet model SICOPOLIS. Two of them are the pre¬ 
viously existing polythermal two-layer (POLY) and cold-ice (COLD) schemes, while the other two are the newly 
implemented conventional one-layer enthalpy (ENTC) and melting-CTS one-layer enthalpy (ENTM) schemes. The 
ENTC scheme goes back to the study by Aschwanden and others ( 2012| l, while the ENTM scheme was introduced 
by Blatter and Greve (20151 for a one-dimensional ice slab. Extension to the full, three-dimensional problem was 
straightforward because, due to the neglect of horizontal diffusive heat fluxes, the computation of the thermodynamic 
fields (temperature and water content or enthalpy) is essentially one-dimensional for each column. Horizontal advec- 
tion merely plays the role of an additional source term for the vertical profiles and does not constitute a problem for 
the implementation of the scheme. 

We used two scenarios from the suite of EISMINT (European Ice Sheet Modeling INiTiative) Phase 2 Simplified 


Geometry Experiments (Payne and others 20001 as our model problems, one without and one with basal sliding. Like 


in the one-dimensional study carried out by Blatter and Greve (20151, the results computed with the POLY scheme 
were assessed to be most reliable and thus chosen as a reference. 

As it was already reported by |Greve| ( 1997b| l, the COLD scheme strongly overpredicts thicknesses of temperate ice 
layers, and thus temperate ice volumes. Choosing a small time step reduces the amount of overprediction to a certain 
extent, but the performance of the COLD scheme (which is not energy-conserving and thus physically incorrect 
anyway) is not satisfactory. 


The studies by Kleiner and others (2015 i and Blatter and Greve (2015 i, both carried out for a Canadian-type 
parallel-sided slab, found that the ENTC scheme, even though it does not enforce the transition conditions at the CTS 
explicitly, still fulfills them for the case of melting conditions (continuity of the enthalpy and the sensible heat flux), 
provided that the discontinuity of the enthalpy diffusivity at the CTS is properly accounted for in the discretization 
of the diffusion term in the enthalpy equation ( [T0) i. This could not be confirmed in the present study; even though 
the ENTC scheme has been implemented in SICOPOLIS in the same way as described by Blatter and Greve| ( |2015| l, 
the scheme fails to produce a continuous temperature gradient across the CTS and overpredicts temperate ice layer 
thicknesses, albeit to a smaller degree than the COLD scheme (always compared to the references provided by the 
POLY scheme). We speculate that this is a consequence of the three-dimensionality of the problem investigated here, 
most likely due to the additional horizontal advection terms in the enthalpy equation. 

In contrast to the ENTC scheme, the continuity of the temperature gradient across the CTS is explicitly enforced 
by the ENTM scheme. This leads to a better match of computed temperate ice layer thicknesses to those computed by 
the POLY scheme; however, some overprediction remains. The price to pay for enforcing a continuous temperature 
gradient across the CTS, while the positioning of the CTS is limited by the grid resolution, is a slight discontinuity of 
the ice temperature itself. 

When comparing the temperature profiles in the ice column above the CTS, the ENTC and ENTM schemes both 
perform well compared to the POLY scheme, while the temperatures produced by the COLD scheme are too high. 
The water content in the temperate ice layer below the CTS usually reaches the threshold value of 0.01 (1%) for the 
POLY, ENTC and ENTM schemes, while the COLD scheme does not compute any water content. 
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For large time steps in conjunction with a water-content-dependent rate factor, the POLY scheme can produce a 
wavy distribution of the temperate-layer thickness. Therefore, some care is required to choose a sufficiently small 
time step to avoid this instability. However, in contrast to the “academic” EISMINT application considered here, real- 
world applications require small time steps for overall numerical stability anyway. At 10 km horizontal resolution, 
Greenland simulations with SICOPOLIS typically require time steps of ~ la, and Antarctica simulations with coupled 
sheet-shelf dynamics even smaller ones. So this instability is most likely not a matter of concern for most practical 
applications with complex topographies. 

The ENTC and ENTM schemes produce some temporal and spatial noise for the volume and thickness of the 
temperate ice. This does not happen for the POLY and COLD schemes; it is related to the one-layer approach in 
which we position the CTS at the uppermost grid point in the temperate layer, whereas the POLY scheme allows, in 
principle, tracking the position of the CTS to arbitrary accuracy. A possible improvement of the one-layer enthalpy 
schemes would be to implement a sub-grid tracking of the CTS position, which would likely allow to fulfill the 
continuity of the temperate and of the temperature gradient across the CTS simultaneously and reduce the noise. 
However, such a sub-grid tracking scheme will inevitably increase the computational cost. 

To sum up, for a polythermal ice sheet model, the one-layer enthalpy schemes ENTC and ENTM are viable, 
easier to implement alternatives to the POLY scheme with its slightly awkward need to handle two different numerical 
domains for cold and temperate ice. ENTM is more precise than ENTC for determining the position of the CTS, and 
thus the volume and thickness of the temperate ice layer. The performance of both schemes is good for the temperature 
profiles in the cold ice column above the CTS and also the water content in the temperate ice layer below the CTS. 
The computing times of all three schemes are comparable. 
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